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ABSTRACT 

Recent proposals have been advanced to apply Imaging Air Cherenkov Telescope 
(IACT) arrays to Stellar Intensity Interferometry. Of particular interest is the possi- 
bility of model-independent image recovery afforded by the good (u, w)-plane coverage 
of these arrays, as well as recent developments in phase retrieval techniques. The capa- 
bilities of these instruments used a s SII receivers have already been explored for simple 
stellar objects (jNunez et alj|2012 ). and here the focus is on reconstructing stellar im- 
ages with non-uniform radiance distributions. We find that hot stars (T > 6000° K) 
containing hot and/or cool localized regions (AT ~ 500° K) as small as ~ 0.1 mas can 
be imaged at short wavelengths (A = 400 nm). 

Key words: instrumentation: interferometers, techniques: high angular resolution, 
techniques: image processing, techniques: interferometric, stars: imaging 



1 INTRODUCTION 

Stellar Intensity Interferometry (SII) has seen a revival 
due to the extraordinary (u, v) covera ge that future air 
Cher enkov telescope arrays will provide (jHolder fc LeBohed 
2006). The angular resolution that can be achieved is 
as fine as 0.06 mas at the longest baselines (1.4 km) and 
the s hortest optical wavelengths (~ 400 nm) ( Nunez et al.l 
The possible improvement in angular resolution 
by an order of magnitude and increased sensitivit y for 
hot (> 6000° K) stellar objects jPravins et alj|20lj ), has 
motivated the exploration of SII capabilities to investigate 
several science topics. These include diameter measure- 
ments, stellar rot ation, gravity darke ning, mass loss and 
mass transfer (see IDravins et al l (|2010h for more details). 

Intensity interferometry consists in measuring the 
squared modulus of the complex mutual degree of co- 
herence between detector pairs from the cross-correlation 
between the light intensity fluctuations. Therefore, only 
the magnitude of the Fourier transform of the stellar 
radiance distribution is accessible. This poses a challenge 
for performing model-independent imaging. However, recent 
developments in phase recovery techniques from simulated 
SII data will make model independent imaging po ssible 
jNunez et al.l 120121 : iHolmes. Nunez fc LeBohed l2010l ). In 
previous work, we have shown that over-all shapes and di- 
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mensions can be reconstructed with sub-percent accuracies 
by using a Cauchy- Riemann phase reconstruction algorithm 
(N uncz et alJ 12012). The use of this phase reconstruction 
algorithm along with iterative post-processing routines also 
allows for further details to be imaged. 

This paper is organized as follows: Section 2 outlines the 
generation of simplified pristine stellar images and the simu- 
lation of intensity-interferometry data. Then phase recovery 
is briefly outlined and image post-processing is presented in 
Section[3] Results are then presented in Section|4]for images 
with increasing degrees of pristine image complexity. 



2 PRISTINE STELLAR IMAGES AND DATA 
SIMULATION 

Pristine images of disk-like stars are first generatecQ. 
These images correspond to black-bodies of a specified 
temperature containing an arbitrary number of "star spots" 
of variable size, temperature, and location at the surface 
of the spherical star in this three dimensional model. The 
simulated stellar surface is then projected onto a plane, so 
that spots located near the edge of the visible half-sphere 
appear more elongated than those located near the center. 



1 The original "pristine" image consists of 2048 X 2048 pixels 
corresponding to ~ 10 mas X 10 mas of angular extension and a 
wavelength of A = 400 nm. 
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Additionally, limb-darkening is included by assuming that 
the stellar atmosphere has a constant opacity. More details 
of the simulated images are presented in section [4] 

From these pristine images, simulated SII array data 
are obtained by computing the Fourier transform of the 
"pristine" image so that the squared modulus of the degree 
of coherence can be found between telescope pairs in 
a large IACT array 0. Gaussian noise is then added in 
the simulations. The signal-to-noise depends on squared 
modulus of the degree of correlation, the area of each 
of the light receivers, the spectral density (number of 
photons per unit area per unit time per unit frequency), 
the quantum efficiency, t he electronic bandw idth, and 
the observation time (see iNufiez et all (|2012h for more 
details.). More detailed simulations which include finite 
detector size, photodetector pulse shape and excess noise 
are currently under development (Rou, Nunez & LeBohec, 
in preparation). 

The simulations used here do not include the change in 
baseline due to earth's rotation, and yet exposure times are 
typically of the order of 10 hrs. These simulations can be re- 
garded as corresponding to the accumulation of a sequence of 
short (< 1 hr) observations taken over several (~ 10) nights. 
It is also possible to save the recorded signal histories in se- 
quences as short as 5 minutes or less, and then this data can 
be correlated in post-processing to avoid smearing of infor- 
mation in the (u, v) plane. The simulated stellar images have 
corresponding temperatures of ~ 6000° K, and the total ex- 
posure time can be significantly reduced if the source has 
a higher temperature. In fact, most stars for which photon 
correlations can be detected within ~ 100 hrs with a large 
IACT array , have temperatures greater than ~ 6000° K (see 
section 2 in lNunez et alj (|2012l ) for details). 



3 CAUCHY-RIEMANN PHASE RECOVERY 
AND POST-PROCESSING ROUTINES 

The Cauchy-Riemann phase recovery algorithm con- 
sists in using the theory of analytic functions to 
relate th e magnitude and the ph ase of the Fourier 
transform (|Holmes fc Belen'kii I 120041 ). That is, in 1- 
dimension, the Cauchy-Riemann equations relate the 
magnitude and the phase differentials along the real and 
imaginary directions in the complex plane. For details 
on the two-di mensional Cauchy -Riemann phase recovery 
algorithm see |Nunez et al.ll2012i ), section 5. The resulting 
reconstructed image with this estimated phase is sometimes 
not ideal, and so is taken as a first guess for the iterative 
algorithms that are now described. 

The Gerchberg-Saxton algorithm 

l|Gerchberg fc Saxton I ll972T ). also known as the error- 
reduction algorithm, is an iterative procedure. Starting 
from a reasonable guess of the image, the algorithm consists 
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Figure 1. Schematic of the Gerchberg-Saxton error reduction 
algorithm. 



in going back and forth between image and Fourier space, 
each time imposing general restrictions in each domain. 
Figure [T] describes the Gerchberg-Saxton algorithm. Start- 
ing from an image Ok, the first step consists in taking 
the Fourier transform to obtain something of the form 
Mh^ k . Next, Fourier constraints can then be applied, 
i.e. the magnitude is replaced by that given by the data, 
and the phase of the Fourier transform is maintained. 
Next, the inverse Fourier transform is calculated and 
constraints can be imposed in image space. The constraints 
in image space can be very general. The image constraint 
that we impose comprises applying a mask, so that only 
pixels within a certain region are allowed to have positive 
non-zero values. For the images presented below, the mask 
is a circular region whose radius is typically found by 
measuring the radius of the first guess obtained from the 
Cauchy-Riemann approach. In all reconstructions where 
the Gerchberg-Saxton is used, we perform 50 iterations, 
and found that more iterations typically do not produce 
significant changes in the reconstruction. 

Another post-processing application that has been uti- 
lized is MIRA (Multi-aperture image reconstruction algorithm) 
jThiebautl 12009). MIRA has become a standard tool for image 
reconstruction in amplitude (Michelson) interferometry. MIRA is 
an iterative procedure which slightly modifies image pixel values 
so as to maximize the agreement with the data. In the image re- 
construction process, additional constraints such as smoothness 
or compactness can be applied simultaneously, but this is some- 
thing that we have not yet experimented with, i.e. the regular- 
ization parameter is set to zero for all reconstructions presented 
here. In the results presented below, the number of iterations is 
set by the default stopping criterion of the optimizer. The MIRA 
software only uses existing data in the (u, v) plane, as opposed to 
using the fit of an analytic function to the data as is done in the 
Cauchy-Riemann and Gerchberg-Saxton routines. This results in 
removing artifacts in the reconstruction th at can be caused by 
the fit of an analytic function to the data. iRenard et al. (2011) 
have recently investigated imaging capabilities with MIRA in the 
context of amplitude interferometry, and the results shown below 
are presented in such a way so that they can be compared with 
amplitude interferometry capabilities. 



2 A preli minary design of the Cherenkov telescop e array (CTA) 
wa s used JCTA Consortium lliuiol; iBernlohr Il2008h . See Figure 2 
in iNufiez et al.ll2012l) for the exact array configuration used in 
these simulations. 



4 RESULTS 

Pristine images are generated with varying complexity. We first 
generate pristine images of stars with limb-darkened atmospheres 
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and investigate the reconstruction capabilities. Then we introduce 
a localized bright or dark feature, and finally increase the number 
of features and explore some of the parameter space, i.e. spot size, 
location, etc. 

4.1 Limb-darkening 

Image reconstruction is actually not necessary for the study of 
limb-darkening, which can be studied with the knowledge of the 
squared degree of coherence only. A direct analysis of the data, 
with no image reconstruction, is likely to yield better results than 
the ones presented below. However, it is instructive to first see 
this effect in reconstructed images before adding stellar features 
to the simulated pristine images. Limb darkening can be approx- 
imately mo deled with a single parameter a as I((fr)/Io = (cos 4>) a 
faestroffeil 1 1997ft . where <f> = refers to light being emitted 
from the center of the stellar disk, and directed radially to the 
observer. The values of a depend on the wavelength and can be 
found by assumi ng hydrostatic equilibrium. At a wavelength of 
400 nm, a « 0.7 l|Hestroffer fc MagnarJll998ft for Sun-type stars, 
and deviations from this value may be indicative of stellar mass 
loss. An example of the reconstruction of a limb darkened star 
with a = 5 is shown in Figure [2] such a large value is chosen 
so that the effect is clearly visible in a two dimensional image 
with linear scale. More realistic values are considered below. 
To obtain Figure [2] a first estimate of the phase was obtained 
from the Cauchy-Riemann algorithm and used to generate a raw 
image. Then the Gerchberg-Saxton post-processing loop (Figure 
[TJ was performed several (50) times. In Figure[3]the ratio of the 
average radius at half maximum R1/2 and the nominal radius of 
Ro, is shown as a function of the limb darkening parameter a. 
Here data were simulated corresponding to stars with apparent 
visual magnitude m v = 3, temperature T = 6000° K, radii 
of Ro = 0.3 mas, lOhrs of observation time and A = 400 nm. 
The ratio R1/2/R0 is less than 1, even in the absense of noise 
since the reconstruction is at best a convolutior[£| with the array 
point-spread- function (PSF). 

From Figure [3] we can see that we are sensitive to changes of 
order unity in the limb-darkening parameter a. Even though this 
analysis is not sensitive to small changes in a, stars experiencing 
high mass loss rates are likely to have high values of a. If we fit 
a uniform disk function to a limb-darkened reconstruction, the 
fit yields a smaller radius. For example, in the case of a = 2.0, 
a uniform disk fit yields an angular radius that is smaller by 7% 
(still larger th an the sub-percent uncertainties found in radius 
measurements jNunez et al,|[2012ft ). A real example is the case 
of the star Deneb, where the difference between the extracted 
uniform disk diameter (Qjjd = 2.40 ± 0.06 mas) and the limb- 
darkened diameter is 0.1 mas l|Aufdenberg et aljliooll . and its 
measured mass loss rate is 10 — 7 M0yr — 1 . 

4.2 Stars with single features 

Stars were simulated as black bodies with a localized feature 
of a higher or lower temperature as described in section [2] In 
the simulated images, the effect of limb darkening is included as 
described in the previous section. Here the full reconstruction 
analysis was used, which consists in first recovering a raw image 
from the Cauchy-Riemann algorithm. The raw image is then used 
as a starting point for several iterations of the Gerchberg-Saxton 

3 A "perfect" reconstruction gives R1/2/R0 = 0.96 for a = 0. 
The extra discrepancy is due to a small hot-spot in the recon- 
struction. R1/2 < Ro since the reconstruction is normalized to 
the highest pixel value. 
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Figure 2. Image reconstruction of a star with limb darkening pa- 
rameter a = 5, apparent visual magnitude m v = 3 and lOhrs of 
observation time. The pristine starting image from which inten- 
sity interferometric data were simulated is shown in the bottom 
left corner with the same contour lines. The Cauchy-Riemann 
phase reconstruction was performed to produce a raw image, and 
then the Gerchberg-Saxton routine was implemented to produce 
the post-processed image shown. 
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Figure 3. For each value of a, SII data were simulated corre- 
sponding to stars with apparent visual magnitude m v = 3, lOhrs 
of observation time and A = 400 nm. For each image reconstruc- 
tion, a 1-dimensional profile is found by calculating the average 
intensity at each radial position. The half height R\/2 (angular 
radius at half maximum) is found. For each value of o, the process 
was repeated 10 times, and error bars are the standard deviation 
of the distribution of H 1 /2(a). The ratio R1/2/R0 < 1 since the 
image reconstruction is at best a convolution with the PSF of the 
array. 

loop (see Figure [TJ, and finally the output of the Gerchberg- 
Saxton algorithm is the starting image for the MIRA algorithm. 
Examples can be seen in Figures [4] and [5] corresponding to the 
post-processed reconstructions of bright stars of m v = 3, 10 hrs 
of observation time and a temperature T = 6000° K. In Figure [4] 
the temperature of the spot is T S p t = 6500° K, and in Figure [5] 
the temperature of the spot is T apo t = 5500° K. 

We can estimate the smallest temperature contrast that can 
be detected by varying the parameters in the model producing 
the pristine image. The performance in terms of temperature 
contrast obviously depends on several variables such as the 
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Figure 4. Reconstructed bright spot. This simulated reconstruc- 
tion corresponds to a star of m v = 3, 10 hrs of observation time, 
T = 6000 K, and spot temperature of T spoi = 6500° K. 
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Figure 5. Reconstructed dark spot. This simulated reconstruc- 
tion corresponds to a star of m v = 3, 10 hrs of observation time, 
T = 6000 K, and spot temperature of T apot = 5500° K. 



size, location and shape of the spot. To quantify the smallest 
detectable spot temperature contrast, we calculate the Mean 
Square Error (MSE] 4 ] comparing the reconstructed image and 
the pristine image convolved with the array PSF. The ratio of 
the angular size of the image and the region over which the 
MSE is calculated (0. 6 mas/1.76 ma s ) is c omparable to the same 
ratio in the work of iRenard et all J201lf> (34 mas/100 mas), so 
that the MS Es presented below a re directly comparable to those 
obtained by IRenard et al.l j201lT> . The MSE is still difficult to 
interpret by itself, so it is compared with the MSE between the 
reconstruction and a simulated star with no spots. By comparing 
these two values we can have an idea of the confidence level for 
reconstructing spots with different temperatures. This compar- 
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Figure 6. The bottom curve and data points correspond to 
the MSE between reconstructed images and their corresponding 
pristine images containing spots of temperature T + AT, where 
T = 6000° K (m„ = 3, 10 hrs of observation time). The top curve 
and data points are the MSE comparing the reconstructed image 
and a uniform disk of the same size as the pristine image. To 
estimate the uncertainties, we performed several reconstructions 
and found the statistical standard deviation of the MSE. 



ison is shown in Figure [6] where the bottom curve corresponds 
to the MSE as a function of spot temperature 6000° K + AT 
between the reconstructed images and the pristine images, and 
the top curve is the MSE comparing the reconstructed images 
and a spotless disk of the same size as the pristine image. A 
total of 26 stars were simulated, and the uncertainty in the MSE 
was estimated by performing several (10) reconstructions for 
one particular case (AT = 500° K). From the Figure it can be 
seen that spots are accurately imaged when AT < — 700° K or 
AT > 200° K approximately. For a black body of spectral density 
B(T), a temperature difference AT < —700° K corresponds to a 
flux ratio B(T + AT) / B(T) < 0.45, and a temperature difference 
AT > 200° K corresponds to flux ratios B(T + AT)/B(T) > 1.2. 
This asymmetry can be partly understood in terms of the 
brightness ratio between black bodies B(T + AT)/ B(T), whose 
rate of change is higher when AT > than when AT < 0. This 
however does not fully account for the asymmetry between cool 
and hot spots. Most of the asymmetry is due to the fact that 
all the simulated stars in Figure [6] have the same integrated 
brightness, and the radiance per solid angle is larger for a star 
containing a bright spot than for an annular region in a star 
containing a dark spot. The same analysis can be performed 
by simulating stars with different integrated brightness, but 
the estimate becomes unnecesarily cumbersome and implies 
knowledge that we would not have access to prior to performing 
an image reconstructiorQ We should also not forget that this is 
an estimate, and in a more precise calculation we would need to 
consider additional variables such as spot size, position, etc. 

To test whether the full chain of algorithms is needed to 
produce Figures [4] and [5] and to investigate algorithm perfor- 
mance, we reconstruct Figs. [4] and [5] with different algorithms 
and combinations of algorithms. Then we calculate the MSE of 
the reconstructions with the pristine image (convolved with the 
PSF), and also compare it to the MSE with a uniform disk. The 
MSE is found for reconstructions using combinations of Cauchy- 



4 For two images Ai j and Bi j of dimension N X N, MSE = 
Min M jiV -2 Y^f^'J ~ B i+k,j+l) 2 ^- Images should be nor- 
mahzed, i.e. y ^ , Ai : j = y t . . Bij = 1. 



5 For example, we would need to have information on the radi- 
ance per solid angle in an annular region in a star containing a 
dark spot. 
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Figure 7. Image reconstructions for different algorithm combi- 
nations. The pristine corresponds to that of Figure [4] 



Riemanrifl (CR), Gerchberg-Saxton (GS), and MIRA. Whe only 
MIRA or Gerchberg-Saxton are used directly (not preceded by 
Cauchy-Riemann), the initial guess of the image is a uniform disk. 
The results are shown in Table[T] The resulting reconstruction for 
each of the algorithm combinations is shown in Figures [7] and [8] 
The single most effective algorithm for these reconstructions is 
MIRA, and the lowest MSE is obtained by using Cauchy Rie- 
mann, followed by Gerchberg-Saxton and MIRA. According to 
table [T] the GS algorithm by itself does better than the CR algo- 
rithm by itself. This is because the pristine images used here have 
a high degree of symmetry. This preferentially aids the GS algo- 
rithm. Also the GS algorithm does better on these cases because 
the support constraint, given by the size of the star, is easy to 
compute. When MIRA or the Gerchberg-Saxton algorithms are 
used directly, the reconstructed image is usually symmetric. The 
Cauchy-Riemann stage provides sensitivity to asymmetries such 
as with the bright or dark spot displayed in Figures [4] and [5] The 
role of Gerchberg-Saxton is more to improve the phase reconstruc- 
tion. MIRA plays the important role of removing artifacts, caused 
for example by the data fitting in the Cauchy-Riemann phase, and 
improving overall definition. Even though non-symmetric images 
can be reconstructed, the final product still is somewhat more 
symmetric than the pristine image for reasons that are still under 
investigation. When the MSE with the pristine image is compared 
to the MSE with a uniform disk, we can again see that the bright 
spot (Figure |4]l is more easily detected than the dark spot (Figure 

EJ. 

A better estimate of the smallest temperature difference that 
can be imaged requires an exhaustive exploration of parameter 
space, but temperature differences of less than a few hundred 
Kelvin do not seem to be possible to image when the same bright- 
ness, temperature, exposure time, angular diameter, spot size and 
spot position as above are used. Results are likely to improve for 
hotter stars than those simulated above since signal-to-noise is 
higher and also the brightness contrast is higher for the same rel- 
ative temperature differences (AT/T). Another question related 
to imaging single spots is that of finding the smallest spot that 
can be reconstructed. In previous work l lNunez et alj|2012l) , we 
show that the smallest possible spot that can be reconstructed 
is given by the PSF of the IACT array used in the simulations, 
namely 0.06 mas. 



6 It only makes sense to use the Cauchy-Riemann algorithm first, 
since this is not an iterative algorithm relying on a first guess. 
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Table 1. Mean square errors for different combinations of recon- 
struction algorithms starting from simulated data corresponding 
to the pristine image of Figs. [4] (top table) and [5] (bottom) . The 
algorithms are CR (Cauchy-Riemann), GS (Gerchberg-Saxton) 
and MIRA. The MSE between the reconstructed image and the 
pristine image (convolved with the PSF of the array) is shown in 
the second column and is compared to the MSE comparing the 
reconstructed image an a uniform disk in the third column. The 
uncertainty in the MSE is AMSE = 4 x 10" 11 . 




Figure 8. Image reconstructions for different algorithm combi- 
nations. The pristine corresponds to that of Figure [5] 



4.3 Multiple features 

As a natural extension to the simulations presented above, 
data were simulated corresponding to stars with two or more 
recognizable features. In Figures [9] and 1101 reconstructions of 
stars containing several hot spots are shown. The brightness 
and exposure time are the same as those used to simulate 
single-spot stars (m v = 3, T = lOhrs). A detailed investigation 
of reconstruction of two-spot stars was not performed, but the 
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Figure 9. Reconstructed star with two hot spots. The pristine 
image has a temperature of 6000° K, and each hot spot has a 
temperature of 6500° K. The simulated data corresponds to m v = 
3 and T = 10 hrs. 

general behavior is similar to that presented in the section 14.21 
The reconstructions improve significantly when the pristine 
image has a higher degree of symmetry, e.g. when both spots lie 
along a line that splits the star in two. This is expected since 
phase reconstruction is not really necessary for centro-symmetric 
images. For this reason, we tested reconstructions with non- 
symmetric pristine images. Even though the shape of the spots 
is usually not well reconstructed, the approximate position and 
size are reasonably accurate. 

The reconstruction and identification of features degrades as 
the number of features in the pristine image is increased. A com- 
mon characteristic of reconstructing stars with several features, 
is that the larger features in the pristine image are better recon- 
structed. This is more so in the case of stars containing darker 
regions. Nevertheless, information of positions, sizes and relative 
brightness of star spots can still be extracted. In Figure ITOl a re- 
construction of a star containing three hot spots of different sizes 
and relative brightness is shown. This simulated reconstruction 
corresponds to the same brightness and exposure parameters as 
those of Figure [9] 



5 CONCLUSIONS 

The capabilities of future IACT arrays for reconstructing complex 
images were demonstrated via simulations. Stars were simulated 
as black bodies with temperatures of 6000° K with localized hot 
or cool regions, then data were simulated as would be obtained 
with an IACT array of ~ 100 telescopes separated by up to 
1.4 Km. A first raw image obtained with the Cauchy-Riemann 
phase reconstruction algorithm was performed, followed by 
post-processing with the Gerchberg-Saxton and MIRA routines. 
Post-processing significantly improves the image reconstruction, 
but the post-processing routines by themselves are usually not 
sufficient for performing reconstructions, especially when the 
pristine image is not centro-symmetric. Features with sizes of 
~ 0.05 mas and temperature differences of a few hundred °K were 
reconstructed. Some statistical tests were performed to quantify 
the confidence with which these features were reconstructed. The 
resulting mean square errors of the simulated reconstructions are 
comparable to simulated results obtained by IjRenard et al]|2011fl 
in the context of amplitude interferometry, but in a different 
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Figure 10. Reconstructed star with three hot spots. The pristine 
image has a temperature of 6000° K, and the spots have tempera- 
tures of 6500° K (top right spot and left spot) and 6800° K (lower 
spot), The simulated data corresponds to m v = 3 and T = 10 hrs. 



angular resolution regime. 

Results suggest that a variety of physical phenomena can be 
imaged. For example, a hot star that is losing mass radiatively 
can have localized bright and dark regio ns corresponding to a 
higher or lower local mass loss rate jLepine fc Mo ffat 2008). The 
origin of these localized regions is unknown, but is likely related 
to magnetic activity. An order of magnitude estimat^ shows that 
stars with a local radiative flux that is roughly twice as high as 
the rest of the star, such as those shown in the Figures above, 
undergo a localized mass loss rate that is three times as high 
as the rest of the star. A related effect that is within imaging 
reach is the von-Zeipel effect, in which the local radiative flux is 
proportional to the surface gravity, so that rotationally distorted 
stars have brighter poles. The analysis of featured oblate star 
observations would be very similar to the one done here. From 
previous results obtained in (Nunez et al. 2012), oblateness ratios 
(ratio of semi-major and semi-minor axis) of 10% — 20% are within 
reach. These in turn correspond to flux ratios of the order of 
20% — 40%, which have also been simulated and reconstructed 
herein. A more comprehensive study on the imagability of these 
physical phenomena will be the subject of future research. It will 
be interesting to see how well radiative-transfer models can be 
tested with the proposed approach. 
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